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Active colloids exhibit persistent motion, which can lead to motility-induced phase separation 
(MIPS). However, there currently exists no microscopic theory to account for this phenomenon. 
We report a first-principles theory, free of fit parameters, for active spherical colloids, which shows 
explicitly how an effective many-body interaction potential is generated by activity and how this 
can rationalize MIPS. For a passively repulsive system the theory predicts phase separation and pair 
correlations in quantitative agreement with simulation. For an attractive system the theory shows 
that phase separation becomes suppressed by moderate activity, consistent with recent experiments 
and simulations, and suggests a mechanism for reentrant cluster formation at high activity. 
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Active colloidal particles in suspension are currently 
the subject of considerable attention, due largely to their 
ability to model self-organisation phenomena in biolog¬ 
ical systems, but also as a new branch of fundamental 
research in nonequilibrium statistical mechanics: assem¬ 
blies of active colloids are intrinsically out-of-equilibrium 
systems. In contrast to their passive counterparts, ac¬ 
tive colloids undergo both solvent induced Brownian mo¬ 
tion and a self-propulsion which requires a continual con¬ 
sumption of energy from the local environment. Several 
idealized experimental model systems have been devel¬ 
oped, such as catalytic Janus particles m. colloids with 
artificial flagella [1] and light activated particles [S]. The 
understanding of active systems has been further aided 
by the development of simple theoretical models, which 
aim to capture the essential physical mechanisms and 
which have been used to study e.g. bacteria, cells or 
filaments in the cytoskeleton [6H9]. 

Active particles are characterised by a persistent mo¬ 
tion, which can lead to ‘self trapping’ dynamics and a rich 
variety of related collective phenomena [MIQ]. Even the 
simplest models of active spherical particles with purely 
repulsive interactions can display the phenomenon of 
motility-induced phase separation (MIPS) (Tn]. In many 
respects, MIPS resembles the equilibrium phase sepa¬ 
ration familiar from passive systems with an attractive 
component to the interaction potential (e.g. the Lennard- 
Jones potential) [TTVH5) . This apparent similarity has 
motivated several recent attempts to map an assembly 
of active particles onto a passive equilibrium system, in¬ 
teracting via an effective attraction (usually taken to be 
a very short range sticky-sphere potential [MIIIT]). De¬ 
spite the intuitive appeal of mapping to an equilibrium 
system, there exists no systematic theoretical approach 
capable of predicting an effective equilibrium potential 
directly from the bare interactions. 

Our current understanding of MIPS has largely been 
gained through either simulation [Mil [IH] or phe¬ 
nomenological theory [ini ni m [IS]. The phenomeno¬ 
logical theory is based on an equation for the coarse¬ 
grained density, featuring a local speed and a local ori¬ 
entational relaxation time. Although the precise rela¬ 


tionship between these one-body fields and the interpar¬ 
ticle interaction potential remains to be clarified, some 
progress in this direction has been made |20j . On a more 
microscopic level, it has recently been shown that a gen¬ 
eral system of active particles does not have an equation 
of state [21] , due to the influence of the confining bound¬ 
aries, however, one can be recovered for the special case 
of active Brownian spheres mm- 

Here we report a first-principles theory for systems of 
active Brownian spheres, which demonstrates explicitly 
how an effective many-body interaction potential is in¬ 
duced by activity. An appealing feature of this approach 
is that intuition gained from equilibrium can be used to 
understand the steady-state properties of active systems. 
The required input quantities are the passive (‘bare’) 
interaction potential, the rotational diffusion coefficient 
and the particle propulsion speed. The theory gener¬ 
ates as output the static correlation functions and phase 
behaviour of the active system. For a repulsive bare in¬ 
teraction, activity generates an attractive effective pair 
potential, thus providing an intuitive explanation for the 
MIPS observed in simulations [niiiiiia. For an attrac¬ 
tive bare potential, we find that increasing activity first 
reduces the effective attraction, consistent with the ex¬ 
periments of Schwarz-Linek et al. [16j . before leading at 
higher activity to the development of a repulsive poten¬ 
tial barrier. We speculate that this barrier may be related 
to the reentrant phase behaviour observed in simulation 
by Redner et al. [14j . 

The paper will be structured as follows: In (jljwe spec¬ 
ify the microscopic dynamics and describe how to elimi¬ 
nate orientational degrees of freedom. From the resulting 
coarse-grained, non-Markovian Langevin equation we de¬ 
rive a Fokker-Planck equation for the positional degrees 
of freedom, from which we identify an effective pair po¬ 
tential. In (jnj we employ the effective pair potential in 
an equilibrium integral equation theory and investigate 
the structure and phase behaviour of both repulsive and 
attractive bare potentials. In the former case we pre¬ 
dict MIPS, whereas in the latter case phase separation 
is suppressed by activity. Finally, in §III| we discuss our 
findings and provide an outlook for future research. 
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I. THEORY 


B. Fokker-Planck equation 


A. Microscopic dynamics 

We consider a three dimensional system of N active, in¬ 
teracting, spherical Brownian particles with spatial coor¬ 
dinate ri and orientation specified by an embedded unit 
vector Pj. Each particle experiences a self propulsion of 
speed Vo in its direction of orientation. Omitting hydro- 
dynamic interactions the particle motion can be modelled 
by the overdamped Langevin equations 

r^ = voP,+'f~^F,+ ( 1 ) 
P^ = ViX P^, ( 2 ) 

where 7 is the friction coefficient and the force on par¬ 
ticle i is generated from the total potential energy ac¬ 
cording to Fi = —\/iUN- The stochastic vectors ^^{t) 
and Piit) are Gaussian distributed with zero mean and 
have time correlations = 2Dtl6ij6{t—t') and 

= 2Drl6ijS{t — t'), where Dt and Dr are the 
translational and rotational diffusion coefficients. 

Equations Q and ([^ are convenient for simulation, 
but are perhaps not the most suitable starting point for 
developing a first-principles microscopic theory. For a 
homogeneous system, averaging over the angular degrees 
of freedom generates a coarse-grained equation [ 12 ] 

+ (3) 

where Xii't) is a Markov process with zero mean and 
where the time correlation function is given by 

= (4) 

The average in (|4) is over both noise and initial orien¬ 
tation. The distribution of Xi{t) is Gaussian to a good 
approximation. This point and further technical details 
of the coarse graining are discussed in Appendix]^ Equa¬ 
tion ^ provides a mean-field level of description, which 
deviates from the exact equations Q and ([^ by neglect¬ 
ing the coupling of fluctuations in orientation and posi¬ 
tional degrees of freedom. 

The Langevin equation (§ describes a non-Markovian 
process, which approximates the stochastic time evolu¬ 
tion of the positional degrees of freedom. The persistent 
motion of active particles is here encoded by the expo¬ 
nential decay of the time correlation with persistence 
time Tp = {2Dr)~^. For small Tp the time correlation be¬ 
comes {Xi{t)Xj{t')) = 2DalSij6{t — t') and the dynamics 
reduce to that of an equilibrium system with diffusion 
coefficient Dt + D^, where Da = VQ/{6Dr). This limit 
is realized when Tp is shorter than the mean free time 
between collisions, i.e. in a dilute suspension. To treat 
finite densities requires an approach which deals with per¬ 
sistent trajectories. With this aim, we adopt (§ as the 
starting point for contracting a closed theory. 


A stochastic process driven by colored noise, such as 
that described by equation is always non-Markovian. 
Consequently, it is not possible to derive an exact Fokker- 
Planck equation for the time evolution of the proba¬ 
bility distribution |24] . Nevertheless, an approximate 
Fokker-Planck description capable of making accurate 
predictions can usually be found. The approximate 
Fokker-Planck equation implicitly defines a Markov pro¬ 
cess which best approximates the process of physical in¬ 
terest (although precisely what constitutes the ‘best’ ap¬ 
proximation remains a matter of debate). From the ex¬ 
tensive literature on this subject (see [SlUSH] and refer¬ 
ences therein) has emerged a powerful method due to 
Fox |271128j , in which a perturbative expansion in pow¬ 
ers of correlation time is partially resummed using func¬ 
tional calculus. The resulting Fokker-Planck equation is 
most accurate for short correlation times (‘off white’ noise 
[25] 1 and for one-dimensional models makes predictions 
in good agreement with simulation data |26j . 

We now consider applying the method of Fox Ell [28] 
to equation This approach consists of first formu¬ 
lating the configurational probability distribution as a 
path (functional) integral and then making a time-local, 
Markovian approximation to this quantity. Technical de¬ 
tails of the method are given in Appendix]^ Fox’s ap¬ 
proach was originally developed to treat one-dimensional 
problems [m Eg, however the generalization to three 
dimensions is quite straightforward. This enables us to 
directly obtain the following Fokker-Planck equation 

N 

at4'(r^t) =-^ V, • J,(r^t), (5) 


where 4' (r^ t) is the configurational probability distri¬ 
bution. Within the generalized Fox approximation the 
many-body current is given by 


J,(r^,t) = -A(r^) V,-/3Ff(r^) 'k(r^,t), (6) 


where /3 =(/cbT) The diffusion coefficient is given by 


D,{r^)=Dt + Da 


\ ^ l-TV,-/3F,(r^)y'’ 


(7) 


where we have defined a dimensionless persistence time, 
T = TpDt/cP. The effective force is given by 

Ff{r^) - fcBTV.I?.(r^)), (8) 


where Di{r^) = Diir^)/Dt is a dimensionless diffusion 
coefficient. Either in the absence of interactions or in 
limit of large Dr the diffusivity Q reduces to Dt -I- Da 
and the effective force becomes DtFi{r^)/{Dt + Da). In 
this diffusion limit the system behaves as an equilibrium 
system at effective temperature Tes = T{1 + Da/Dt). 
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For weakly persistent motion, r—>-0, equations ([^ to 
^ become exact and the theory provides the leading or¬ 
der correction to the diffusion approximation. However, 
the Fox approximation goes beyond this by including con¬ 
tributions to all orders in r. Indeed, detailed studies of 
one-dimensional systems have demonstrated good results 
over a large range of r values |55] . The only caveat is that 
the condition 1 —TVi-/3Fi>0 must be satisfied [IZIIIH]. 
The range of accessible r values thus depends upon the 
specific form of the bare interaction potential. 

Within our stochastic calculus approach, the effec¬ 
tive many-body force (|^ emerges in a natural way from 
the coarse grained Langevin equation ([^. The more 
standard route (adopted in all attempts made so far 
[101 US]) to approach this problem is to derive from the 
Markovian equations Q the exact Fokker-Planck equa¬ 
tion for the joint distribution of positions and orienta¬ 
tions, P(r^p^t). However, coarse graining strategies 
based on integration of P over orientations generate in¬ 
tractable integral terms. By starting from ^ we are able 
to circumvent these difficulties. As we shall demonstrate 
below, our effective force accounts for several important 
collective phenomena in active systems. 


C. Effective pair potential 


In the low density limit we need only consider isolated 
pairs of particles. In this limit reduces to an equation 
of motion for the radial distribution function, g{r,t) = 
'I'(r, t)jpl, where ph is the bulk density. This equation of 
motion, the pair Smolochowski equation, is given by 


(9t5(r,t) = -V • j(r,t), (9) 


where r = \ri 2 \ is the particle separation and V = 
The pair current is given by 


j(r,t) = -2D{r)g{r,t) Vlng(r,t) - j3F°^{r) 
where the radial diffusivity 

TV'^/3u(r) 


Dir) =Dt + Da[l- 


I -I- rV^/Suir) 


( 10 ) 


( 11 ) 


interpolates between the value Dt at small separations, 
where M(r) is strongly repulsive, and Dt + Da at large 
separations. The effective interparticle force is given by 

F^^ir) = :^^(F{r)-kBTVV{r)), (12) 

where the bare force is related to the pair potential by 
Fir) = —Vuir). The symmetry of the two-body prob¬ 
lem can be exploited to calculate from (12) an effective 
interaction potential 




(13) 


where P(r) = |F'(r)|. We have thus identified an effective 
interaction pair potential, which requires as input the 
bare potential and the activity parameters r and Da- 



FIG. 1: Activity induces effective attraction. Passive 
potential /3u(r) = r~^^. (a) Increasing Pe (in steps of 8) 

from 0 to 40 generates an effective interparticle attraction. 
Points indicate the potential minima, (b) Radial distribution 
function, gir), from simulation (points) and theory (lines) 
for pb = 0.5 and Pe = 0 to 20 (in steps of 4). Curves are 
shifted vertically for clarity, (c) As in (b), but focusing on 
larger separations for Pe = 4 (squares), 12 (circles) and 20 
(diamonds). Inset: position of the first peak in gir) as a 
function of Pe. (d) Spinodals for r = 0.045 (dot-dashed) to 
0.065 (long dashed) in steps of 0.005. 


II. RESULTS 


A. Motility-induced phase separation (MIPS) 


To illustrate how activity can generate an effective at¬ 
traction in a passively repulsive system we consider the 
non-specific potential j3uir) = In Fig{^ we show 

the evolution of the effective potential for fixed r 
as a function of the dimensionless velocity Pe = Vod/ Dt. 
For Pe> 10 the effective potential develops an attractive 
tail. As Pe is increased the potential well deepens, the 
minimum moves to smaller separations and the radius of 
the soft repulsive core decreases. These trends are con¬ 
sistent with the intuitive picture that persistent motion 
drives soft particles into one another (the soft core ra¬ 
dius reduces) and that they remain dynamically coupled 
(‘trapped’) for longer than in the corresponding passive 
system. Within our equilibrium picture the trapping is 
accounted for by the effective attraction. 


For systems at finite density the pair potential (13) is 


an approximation, because three- and higher-body inter¬ 
actions will play a role (see equation ([^). However, for 
simplicity we henceforth employ the pair potential (13) 


for all calculations, as we anticipate that this will pro¬ 
vide the dominant contribution. Although corrections to 
this assumption can be made, they obscure the physical 
picture and come at the expense of a more complicated 
theory. The validity of the pair potential approximation 
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is very well described. Further comparison for other pa¬ 
rameter values (not shown) suggests that (131 combined 
with the SMSA theory provides an accurate account of 
the asymptotic decay of pair correlations. 

In Fig[^ we show the spinodal lines mapping the 
locus of points for which the static structure factor, 
S{q) = (1 — Pbc{q))~^, diverges at vanishing wavevec- 
tor. Simulations have shown that MIPS is consistent 
with a spinodal instability m- As T is decreased the 
critical point moves to higher values of Pe and to slightly 
higher densities. When compared with the spinodal of a 
standard Lennard-Jones system (e.g. the black curve in 
Figj^) the critical points in Fig{^ lie at rather higher 
values of pb- This suggests that typical coexisting liquid 
densities for MIPS will be larger than those found in equi¬ 
librium phase separated systems, as has been observed in 
simulation mill]. 


FIG. 2: Prom active suppression of phase separation 
to a cluster phase. Passive potential j3u{r)=As{r~^^—r~^) 
with £ = 1.4. (a) Effective potential ( |13[ ) for t — 0.025 and 
Pe = 0 (black) 4, 8,12 (red broken lines) aiad 20, 28, 36 (blue). 
Inset: zoom of the repulsive peak for Pe = 20,28,36. (b) 
Spinodals for Pe = 0, 4, 8,12. Increasing Pe increases £crit, the 
critical value of e. (c) Scrit as a function of Pe (black full line) 
and the locus of points for which the repulsive peak of 
takes the value 0.1 (red dashed) and 0.05 (blue dot-dashed). 
Open (closed) circles indicate points where BD simulation 
find a mixed (phase-separated) state (see hgure]^. Arrow 
indicates path taken in (a), (d) Theory (lines) and simulation 
(symbols) data (shifted for clarity) for g{r) at £ = 0.5, pb = 
0.3 for Pe = 0 (black), 12 (blue) and 20 (green). Dotted lines 
indicate peak positions. 


is justified a posteriori by the comparison with simulation 
for the finite density pair correlations. 

In Fig{^ we show the steady-state (isotropic) radial 
distribution function for /9b = 0.5 for various values of Pe. 
We employ the effective pair potential together with 
liquid state integral equation theory and compare theo¬ 
retical predictions with direct Brownian dynamics sim¬ 
ulation of equations Q and ([^. The integral equation 
theory we employ is the soft mean-spherical approxima¬ 
tion (SMSA) proposed by Madden and Rice m. This 
approximate closure of the Ornstein-Zernike equation is 
known to provide reliable results for the pair stucture of 
Lennard-Jones type potentials. Given the form of the 
effective pair potential shown in Figj^the SMSA would 
seem to be a reasonable choice of closure. Details of the 
integral equation theory and the simulation procedure 
are given in Appendices [C] and [P] respectively. 

We find that as Pe is increased the main peak of g{r) 
grows in height and shifts to smaller separations (see inset 
to Fig[^), reflecting the changes in the effective poten¬ 
tial. In the main panel of Fig[^ we focus on the second 
and third peaks. The quantitative accuracy of the the¬ 
ory in describing the decay of g{r) is quite striking, in 
particular the phase shift induced by increasing activity 


B. Suppression of phase separation 

We next consider the influence of activity on a 
Lennard-Jones system, I3u{r) = A£{r~^^ — r“®). For a 
phase-separated passive system, recent experiments and 
simulations have demonstrated that increasing Pe first 
suppresses the phase separation m and then leads at 
higher Pe to a reentrant MIPS [14]. Schwarz-Linek et al. 
have argued that the suppression of phase separation at 
lower to intermediate Pe occurs in their system because 
particle pairs bound by the attractive (depletion) poten¬ 
tial begin to actively escape the potential well, and that 
this can be mimicked using an effective potential less at¬ 
tractive and shorter-ranged than the bare potential ng. 

To investigate these phenomena we set £=1.4, which 
ensures a phase separated passive state [31] , and consider 
the evolution of the effective potential as a function of Pe. 
In Fig[^ we show that as Pe is increased from zero to the 
value 18 both the depth and range of the effective poten¬ 
tial reduce significantly, consistent with the expectation 
of Schwarz-Linek et al. [ig. Spinodals within this range 
of Pe values, identifying where the static structure fac¬ 
tor diverges at zero wavevector, are shown in Fig[^. As 
Pe is increased the critical point moves to higher values 
of £ (cf. Figs.l and 3 in Ref.[Tg). A passively phase- 
separated system will thus revert to a single phase upon 
increasing the activity. To examine this behaviour in 
more detail we show in Fig[^ the trajectory of the criti¬ 
cal point in the (Pe, e) plane. Above the line there bulk 
densities for which phase separation occurs. 

In order to test the predicted trajectory of the critical 
point we have performed Brownian dynamics simulations 
at a bulk density /96 = 0.4, which lies close to the critical 
density m, for various values of £ and the Pe values 5, 8 
and 10. Visual inspection of the simulation snapshots re¬ 
veals the existence of voids in the particle configurations 
corresponding to a phase separated state (see Figsj^ and 
for snapshots). This visual impression can be made 
more quantitative by calculating the radial distribution 
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FIG. 3: Simulated phase separation, (a) Snapshot of a 
mixed system at t/rs = 40, Pe = 8, e/iksT) = 1.5. (b) 
Snapshot of a phase separating system at the same time and 
Pe = 8, £ = 2.5. (c) The radial distribution function, g{r), for 
£ = 1.5 (black curve), 2 (blue curve), and 2.5 (green curve), 
(d) As in (c) but focusing on larger distances. From the 
snapshots together with the long range behaviour of the g{r) 
we can distinguish between a mixed and a phase separated 
system. The slow decay to the asymptotic value of unity, as 
shown in (d), indicates phase separation. 


function. Phase separating states generate a very charac¬ 
teristic slow decay of g{r) (Figsj^ and|^) which provides 
a useful indicator. The open circles in Figj^ represent 
mixed states, whereas closed circles indicate phase sepa¬ 
rated statepoints. The phase boundary predicted by the 
theory is highly consistent with the simulation data. 


C. Cluster phase 

Returning to Figj^, we find that for Pe> 18 the effec¬ 
tive potential develops a repulsive barrier, which grows 
in height (see inset) with increasing Pe, while the poten¬ 
tial minimum becomes deeper. It is well-known that po¬ 
tentials with a short-ranged attraction and long-ranged 
repulsion (SALR potentials) exhibit unusual equilibrium 
phase behaviour, including clustering and microphase- 
separation [321133] ■ Although the attractive component 
of the potential may favour phase separation, the long 
range repulsion destabilises distinct liquid and gas phases 
and causes them to break up into droplets or clusters. 
This represents a non-spinodal type of phase transition, 
characterized by a divergence in the structure factor at fi¬ 
nite wavevector. The appearance of a repulsive barrier in 
the effective potential suggests that a similar mechanism 
may be at work in passively attractive systems subject 


to high Pe activity. 

In Figj^ we show the locus of points where the ef¬ 
fective potential peak height attains a given value (we 
choose 0.05 and 0.1 for illustration). When these ‘isore¬ 
pulsion curves’ are viewed together with the critical point 
trajectory the resulting phase diagram is very similar 
to that obtained by Redner et al. in their simulation 
study of two-dimensional active Lennard-Jones particles 
(cf. figure 1 in Ref.[T3]). However, a detailed study of 
the connection between the potential barrier and high Pe 
clustering goes beyond the scope of the present work. 


III. DISCUSSION 

In summary, we have shown that systems of active 
spherical Brownian particles can be mapped onto an 
equilibrium system interacting via an effective, activity- 
dependent many-body potential. The only required in¬ 
puts are the bare potential, thermodynamic statepoint 
and the parameters specifying the state of activity. Our 
theory captures the phenomenon of MIPS in repulsive 
systems and provides first-principles predictions for the 
activity dependence of the pair correlations, in very good 
agreement with Brownian dynamics simulation. As far 
as we are aware there exists no other approach capa¬ 
ble of predicting from the microscopic interactions the 
pair correlation functions of an active system. Further 
insight into the steady state particle distribution could 
in principle be obtained by investigating the three-body 
correlations. These could be obtained by employing the 
effective potential in a higher-order liquid state integral 
equation theory (see e.g. |34j and references therein). 

For passively attractive systems the theory rationalizes 
the experimental finding |16] that increasing activity can 
suppress passive phase separation. We find that as Pe 
is increased from zero to intermediate values the mini¬ 
mum of the effective potential becomes less deep, thus 
weakening the cohesion of the liquid phase. To the best 
of our knowledge there is currently no alternative theo¬ 
retical explanation of phase-transition suppression in ac¬ 
tive suspensions. It is an appealing aspect of our theory 
that the suppression of passive phase separation follows 
naturally from the same approach which yields activity- 
induced attraction for repulsive potentials. For high val¬ 
ues of Pe the appearance of a repulsive barrier in the 
effective potential suggests that the reentrant phase sep¬ 
aration observed in simulations m may be interpreted 
using concepts of equilibrium clustering in SALR poten¬ 
tial systems. This will be a subject of future detailed 
investigations. It is known that care must be exercised 
when analysing SALR potentials, as traditional liquid 
state theories can prove misleading [55] . 

A key step in our development is the Fox approxima¬ 
tion [27] , which yields an effective Markovian description 
of the coarse-grained equation Making a Markovian 
approximation automatically imposes an effective equi¬ 
librium, however, we are aware that there exist certain 
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situations for which this breaks down [iniiiii- Establish¬ 
ing more clearly the range of validity of our approach, as 
well as its possible extensions, will be the subject of ongo¬ 
ing study. However, it is already clear that going beyond 
the Markovian approximation will be very challenging. 
Indeed, such a step may not even be desirable. Any kind 
of non-Markovian description would lead inevitably to a 
loss of the effective equilibrium picture and the phys¬ 
ical intuition associated with it. It thus seems likely 
that practical improvements to the present approach will 
retain the Markovian description while seeking to opti¬ 
mize, or improve upon, the Fox approximation for certain 
classes of bare potential. Very recently, Maggi et al. have 
employed an alternative approach to treating stochastic 
processes driven by Ornstein-Uhlenbeck noise [3S]. A 
comparison of their approach with the Fox method em¬ 
ployed here would be very interesting. 

With a view to further applications of our approach, 
we note that there has recently been considerable inter¬ 
est in active suspensions at very high densities [55H55] . 
In particular, it has been found using computer simula¬ 
tions that activity has a strong influence on the location 
of the hard-sphere glass transition, dynamic correlation 
functions, such as the intermediate scattering function, 
and static pair correlations m- Within our effective 
equilibrium framework, increasing the activity of a pas¬ 
sively repulsive system generates an effective attraction. 
We can therefore anticipate that for volume fractions just 
above the glass transition it will be possible to observe a 
reentrant glass transition, namely a melting of the glass 
followed by revitrification, as a function of increasing Pe. 
Moreover, the nontrivial evolution of the effective po¬ 
tential as a function of Pe for attractive bare potentials 
(cf. Figj^) suggests these systems will present a rich vari¬ 
ety of glassy states. Work along these lines is in progress. 

Finally, we mention that a natural generalization of the 
present theory is to treat spatially inhomogeneous sys¬ 
tems in external fields. Recent microscopic studies of ac¬ 
tive particles under confinement (e.g. in a harmonic trap 
[29| ) have provided considerable insight, however none 
of the existing approaches have considered effective in¬ 
terparticle interactions. Inhomogeneous generalization of 
the present theory enables the interaction between MIPS 
and external fields to be investigated on the microscopic 
level. Our preliminary investigations reveal, for exam¬ 
ple, activity-induced wetting at a planar substrate and 
capillary-condensation under confinement. This will be 
presented in a future publication. 
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Appendix A: Coarse-grained Langevin eqnation 


Equation ([^ describes the orientational diffusion of 
an active particle. The corresponding conditional prob¬ 
ability distribution function T(p, 11 pg, to), where t > to, 
obeys a Fokker-Planck equation which can be obtained 
by usual techniques [40], 

^^iP,t\Po,to) = DrR^T{p,t\po,to), (AI) 


where R = {p x Vp) is the intrinsic angular momentum 
differential operator. Eq.(AI) describes nothing but a 


diffusion process on the unit sphere. This problem is well- 
known when studying, e.g., dielectric relaxation in polar 


liquids [4TVl44j . In spherical coordinates, (AI) becomes 




1 d f . a 

sind^dl di} 


I 02 
sin^ d dif'^ _ 


T(H, 11 Ho, to), 

(A2) 


where we have defined H = (r?, ip). 

Assuming that T and its derivatives are continuous on 
the sphere [45] . we expand the probability distribution 
function T in spherical harmonics 


OO I 

T(H, 1 1 Hg, to) = EE Aim{t\^0,to)Yim{^) , (AS) 

l—O rn— — l 


where Yi^ are the spherical harmonics and are co¬ 
efficients encoding the initial condition. We also recall 
that spherical harmonics are eigenvectors of the operator 
R^ (in spherical coordinates), namely that 

R^Yim = -l{l + . (A4) 


Inserting (AS) in 


and using (|A4|) we obtain 


I Ho, to)Yimi^) = 

l,m 


-DrY, I no,to)Yi^in). 

l,m 


(AS) 


Multiplying both sides of (AS) by (H), integrating 
over solid angle and using the orthogonality property, 
I d^Y{:^,{Vl)Yim(yt) = yields 

'^Aim{t I Ho, to) = —Drl{l -f l)Aim(t \ Ho, to) , (A6) 

at 

which has the solution 

Aim{t I Ho, to) = e-^’-'('+i)(‘-‘Aa,^(Ho), (A7) 


where the aim are a new set of coefficients. The proba¬ 
bility distribution is thus given by 

T(H, 11 Ho, to) = y] e-^^'('+i)(*-‘«)ai^(Ho)Y^(H). 

l,m 

(A8) 
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The initial condition, 

T(r2, to I ^o) = ~ ^o) j 


(A9) 


together with the completeness relation of the spherical 
harmonics, 

OO I 

— r^o) — ^ 5] Yi^in)Yr^{no), (Aio) 

/—0 m— — l 

allows the missing coefficients to be identified, 

aUn^)=Y*^{n^). (All) 

The conditional probability distribution is now fully de¬ 
termined as 

OO I 

l—O m— — l 

(A12) 

As t —>■ OO only the terms with I — 0 survive. The steady- 
state distribution function is thus given by 

Teqi^)= lim T(0,t I = (47r)“^ (A13) 

>-oo 


The conditional and equilibrium distributions, (A12) 


and (A13), respectively, can be used to coarse-grain the 
exact Langevin equations Q and ([^. The approach 
taken is to consider the orientation vector p^ (t) attached 
to particle z as a stochastic variable and to provide its 
full statistical characterization. In spherical coordinates 
the orientation vector is given explicitly by 

pit) = iPxit),Pyit),Pzit))'^ 

= {cos (fit) sin-d(t), siny)(t) sind(t), cosi?(t))^, 

(AM) 

where p and d are the azimuthal and polar angles re¬ 
spectively. Using (A13l we have that 


1 


dVt / dUg E/' 




X Y,*oin)YimmYioi^o)Y:^{no) 


where we have expressed the cosine functions in terms 
of spherical harmonics, Yio = \/3/(47r) cosd = Y]*q, and 
used the orthogonality property. Calculations for the xx 
and yy components are performed in the same spirit. We 
thus obtain 


{Pxit)pxito)) = ‘“I = {Pyit)Py{to)), 


(A19) 


whereas off-diagonal components of the correlation ma¬ 
trix are all zero. We can thus conclude that 

{x^it)xAt')) = vl{p,{t)pAt')) = 


It has been shown 


(A20) 

that the probability distribution 
function (A12| can be well approximated by an expres¬ 
sion which generalizes the planar Gaussian function to 
the sphere. The new noise function Xiit) is thus ap¬ 
proximately Gaussian distributed with zero mean and 
exponentially decaying correlations. The coarse-grained 
Langevin equation (§ thus describes a stochastic process 
with additive colored-noise. 


Appendix B: Approximate Fokker-Planck equation 

To derive from (|^ an approximate Fokker-Planck 
equation we apply the functional calculus methods of Fox 
m- We address the one-dimensional case before gen¬ 
eralizing to higher dimension. Consider the stochastic 
differential equation 


x{t) = F{x) + g{x)xit), 


(Bl) 


{Pz{t)) = j dU Te5(U) cos'd = 0, (A15) 

together with analagous results for the x and y compo¬ 
nents 

{Px{t)) = 0 = {Py{t)) . (A16) 

Defining the new stochastic variable by Xiit) = voPiit), 
its first moment is thus given by 

iXiit)) =vo{pM) =0- (A17) 

Calculation of the equilibrium correlation matrix re¬ 
quires th e con ditional probability distribution function 
given by (A12). For example, for the zz component, we 
obtain 

{Pzit)pzito)) =/dD/dDocosdcosdoT(D,t I Dq, to)Tf'eq(Do) 


where F{x) and g{x) may be nonlinear functions in x. If 
g{x) = I, the process is then called additive, otherwise 
it is called multiplicative. The noise function y(t) is by 
definition Gaussian distributed with zero mean. Its sec¬ 
ond moment determines whether it is a white or colored 
noise. As we are interested here in the case of additive 
colored noise we set g{x) = 1. 

In the framework of functional calculus, the Gaussian 
nature of x(t) is expressed by the following probability 
distribution functional, 


P{'X\ = Ne~ ^ -f ‘^'^'x(e)x(d)K(s-s') ^ 


(B2) 




(A18) 


where the function K is the inverse of the x correlation 
function and the normalization constant is expressed by 
a path integral over y, 

The first and second moments of y are given by 

ixit)) = 0, (B4) 

{xit)xis)) = C{t-s). (B5) 








Recalling that the functional derivative may be defined 
according to 


where we have used (BIO) and (B7). Inserting (B13) 




(B6) 


A=0 


we now derive two useful identities. The first concerns 
the functional derivative of the probability distribution 
functional, 


SN 


^P[x] _ _ 

Sx{t) 


■N- 


^-1 f ds'x{s)x(s')K{s-s') 


= ! ds J ds'x{s)x{s')K{s-s') 


Sxit) 

= -P[x] j dsK{t- s)xis), 


(B7) 


where, using (B6) and (B4), it can be easily shown that 


6N/Sx{t) = 0. The second identity demonstrates the 
inverse relation between the functions K and C. The 
second functional derivative of P\x\ yields 

Sxifl^xit) ^ \^Jds'JdsK{t' - s')K{t - s)x{s')x(s) 

(B8) 


where use of (B7| has been made. Using ( |B8[ ) and (B5) 
together with the normalization f D[x\P[xY= 1, leads to 

0 = Id[x]^^^ 

J 6x{f)Sxit) 

= Jds'K{t'-s') jdsK{t-s)C{s-s')-K{t-t')^ 


back into the second term of (B12) and integrating by 
parts gives us 

J D[x]d{y - x{t))P[x]x{t) 

= - Jds'C{t-s')J d[x] (^^S{y - x{t))'^ 

(B14) 

which serves as the exact starting point for Fox’s approx¬ 
imation scheme m- 

In order to progress further we need to calculate 
dx{t)/Sx{s')- Applying the functional derivative with 
respect to x(t') on (Bl) yields a first-order differential 
equation, 

d 6x{t) Sx{t) 

the solution of which is 

dx[s') Jo 

dsF'(rr(s))0(^_g/), 

where 0 is the Heaviside step function, which we define 
here as follows 

1, t > s' 

0(t — s') = ■{ t = s' 

0, t < s'. 


(B9) Using (B16) in (B14) we can rewrite (B12) in an alter- 


which implies that 

JdsK{t-s)C{s-s') = S{t-s'). (BIO) 

The solution to the stochastic process described by 


native form 

= --^[Fiy)Piy,t)] 


(Bll, namely the probability distribution functional for 
is given by the formal expression 


(B17) 

/ which already begins to resemble a Fokker-Planck-type 
P[x]P[x]diy — x(t)) . (Bll) equation. However, because of the non-Markovian nature 

of dsF'{x{s)) appearing in the exponential of (B17), it 
is apparent that a reduction of this term to an expression 
containing P{y, t) is not possible. An approximation is 
required. 

The colored noise of interest here is characterized by 
an exponentially decaying correlation function (A17|. 


Taking the time derivative of (Bll|) yields 


^-P(?/U) = -^[P{y)Piy,t)] 


d_ 

dy 


D[x]S{y-x{t))P[x]xit)- (B 12 ) 


The product P[x\ x(0 appearing in the second term can 
be rewritten in the following way 

P[x] X(J) = P[x] J ds6{t - s)x(s) 


In the literature on non-Markovian processes the time- 
correlation functions are generally notated as follows 


C{t-s) = -e-^ 

T 


(B18) 


with a diffusion coefficient D and a correlation time r. 
In order to retain some coherence with the existing liter¬ 
ature we will here employ the standard notation of (B18) 
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and only use the relation of the parameters in (B18) to 
those of (A20) at the end of the calculation. 

Returning to (B17), we first perform a change of vari¬ 
able, t' = t — s', in the time-integral, 

Jo Jo 

(B19) 

and then expand the time integral over F' in terms of t', 

dsF'{x{s)) « F'{x{t))t' - F"(a;(t))i(t)y. (B20) 

Neglecting the term in (B20) enables the integral in 
(B19) to be evaluated 

[ ds'C[t- ^ 

Jo Jo 

(B21) 


T In 


D 


1 — TF'{x(t)) ’ 


where we used (B18), and the second approximation re¬ 
sults from assuming a sufficiently large t. We can fi¬ 
nally put (B211 back into (B171 to obtain an approximate 
Fokker-Planck equation. 




P{y,t) 

(B22) 


This is Fox’s result for the approximate Fokker-Planck 
equation corresponding to the non-Markovian process 
(Bl|. Eq.(B22) implicity defines a Markovian process, 
which approximates the non-Markovian process of phys¬ 
ical interest. However, the question of whether this rep¬ 
resents the best approximation remains a subject of de¬ 
bate. We note that equation (B22) has also been derived 
by Grigolini et al. |26j using alternative methods which 
do not make any assumptions of a short correlation time. 

The one-dimensional Fokker-Planck equation (B22) 
can be generalized without much difficulty to describe 
a three-dimensional system of N particles. The dynam¬ 
ics of interest is described by the stochastic equation (|^ . 
We now adapt the standard notation used above to that 
employed in the main text, namely P{y,t) —>■ 

T —>■ Tp = l/{2Dr) and D —)■ Ug/S, and recall that 
Da = VQ/iGDr) and = /3Dt for the friction coefficient 
in ([^. Making the appropriate replacements enables us 
to write the three-dimensional generalization of (B22) 

r) ^ 

t) = - ^ V, • A [l3F,{r^) - V,] vE-(r^, t) 


N 




-AV,; 


1 


DoVi-PFijr") 


1 - 




2Dr 


(B23) 


A simple rearrangement of terms in (B231 leads directly 
to equations (©-§ in the main text. 


Appendix C: Integral equation theory 


To calculate the steady-state radial distribution func¬ 
tion, g{r), from the effective pair potential (13) we em¬ 


ploy an equilibrium liquid state integral equation devel¬ 
oped by Madden and Rice [50] . This soft mean-spherical 
approximation (SMSA) exploits the Weeks-Chandler- 
Anderson splitting of the pair potential m into attrac¬ 
tive and repulsive contributions, u{r) = u^epir) -l-Matt(?’), 
where the repulsive part is given by 


-*'rep 


(r) = 


u{r) - u(rmin) 

0 


r < Tmin 

r > r^in 


and the attractive part is given by 

I U{r) r > Tmin 




I u(^niin) 


r <r„ 


(Cl) 


(C2) 


where r-a^in is the position of the potential minimum. The 
total correlation function, h{r) = g{r) — 1, is related to 
the shorter range direct correlation function, c(r), by the 
Ornstein-Zernike equation |35| 


hir) = c[r) + pb J dr'h{\r — r'\)c{r'). 


(C3) 


The SMSA approximation is given by the closure relation 
c(r) = (1 - - /3uatt(?')- (C4) 


For the Lennard-Jones potential the closure relation (C4) 


has been shown to provide results for g{r) which are 
superior to both Percus-Yevick (PY) and Hypernetted 
Chain (HNC) theories [5D|. Moreover, the SMSA the¬ 
ory predicts a true spinodal line in the parameter space, 
namely a locus of points for which the static structure 
factor, S{k) = (1 — pbc{k))~^, diverges at vanishing 
wavevector. This behaviour is a consequence of the as¬ 
sumed asymptotic form of the direct correlation function, 
c(r) ~ —/3uatt (■'’)• Other standard integral equation the¬ 
ories, such as PY and HNC, do not exhibit a complete 
spinodal line, but rather a region within which the theory 
breaks down (‘no solutions region’) [49] . 


Appendix D: Brownian Dynamics Simulations 

In order to benchmark our theoretical predictions we 
perform Brownian dynamics simulations of N particles, 
randomly initialized without overlap. The system is con¬ 
fined to a periodic cubic box, the size of which is deter¬ 
mined by the number density according to = N/pb, 
where L is the side length. The Langevin equations of 
motion, ([^ and ([^, are integrated via a standard Brow¬ 
nian dynamics scheme |50j with a constant time step 
of SI/tb = 10“®. Both the translational and rotational 
noise are Gaussian random variables with a standard de¬ 
viation of (Tt = ( 2 AY )2 and ar = {‘ 2 DrT )2 ^ respectively. 
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For the soft repulsive potential to be considered in this 
work, f3u{r) = we employ N = 2000 particles. The 
potential is truncated and shifted at Vcut/d = 2. To pro¬ 
vide good statistics for the static quantities the simula¬ 
tions are carried out for 10® time steps, sampling every 
1000 steps, which is equivalent to a total run time of 
ttotlTB = 10 and a sampling rate of Ts/tsample = 100. 


For the second system we will consider, the Lennard- 
Jones system, /3u{r) =4e(r“^^—r“®), we simulate a larger 
system of 5000 particles. The integration time of the 
equations of motion is the same as in the repulsive sys¬ 
tem, as is the cut-off radius. In this case, the runtime 
is 10^ and the particle positions are sampled every 10^ 
steps. 
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